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We investigate the dynamics of a neural network where each neuron evolves according to the com- 
bined effects of deterministic integrate-and-fire dynamics and purely inhibitory coupling with K 
randomly-chosen "neighbors". The inhibition reduces the voltage of a given neuron by an amount 
A when one of its neighbors fires. The interplay between the integration and inhibition leads to a 
steady state which is determined by solving the rate equations for the neuronal voltage distribution. 
We also study the evolution of a single neuron and find that the mean lifetime between firing events 
equals f + KA and that the probability that a neuron has not yet fired decays exponentially with 
time. 



PACS numbers: 87.f8.Sn, 87.f9.La, 07.05.Mh. 



I. INTRODUCTION 



Networks of neurons which undergo "spiky" dynamics 
have been thoroughly investigated (see e.g. (TJH] and ref- 
erences therein). Nevertheless, a theory which describes 
the dynamics of randomly interconnected excitatory and 
inhibitory spiking neurons is still lacking. Even a system 
composed exclusively of inhibitory neurons appears 
too complicated for analytical approaches. Part of the 
reason for this is that the dynamics of a single neuron in- 
volves many physiological features across a wide range of 
time scales that arc difficult to incorporate into an ana- 
lytical theory. Our goal in this work is to describe some of 
the dynamical features of a purely inhibitory neural net- 
work within the framework of a minimalist model. While 
we sacrifice realism by this approach, our model is ana- 
lytically tractable. This feature offers the possibility that 
more realistic networks may be treated by natural exten- 
sions of our general framework. 

We specifically investigate an integrate-and-fire neural 
network, in which the integration step is purely linear 
in time, and in which there exists only inhibitory and 
instantaneous coupling between interacting neurons. We 
thus ignore potentially important features such as voltage 
leakage during the integration, as well as heterogeneity 
in the external drive and in the network couplings. How- 
ever, we do not assume all-to-all coupling, an unrealistic 
construction which is often invoked as a simplifying as- 
sumption. Instead, the (average) number of "neighbors" 
is a basic parameter of our model. Analytically, we con- 
sider annealed coupling, where the neighbors of a neuron 
are reassigned after each neuron firing event. However, 
our simulations indicate that the model with quenched 
coupling, where the neighbors of each neuron are fixed 
for all time, gives nearly identical results. 

Each integrate-and-fire neuron is represented by a sin- 
gle variable - the polarization level, or voltage V. Our 
model has two fundamental ingredients: (i) the dynam- 
ics of individual neurons, and (ii) the interaction between 
them. For the former, we employ deterministic integrate- 



and-fire dynamics, in which the voltage on a single neu- 
ron increases linearly with time until a specified thresh- 
old is reached B . At this point the neuron suddenly fires 
by emitting an action potential and the neuronal voltage 
quickly returns to a reference level (Fig. |l|) . For concrete- 
ness and without loss of generality, we assume the rate 
of voltage increase and the threshold voltage are both 
equal to 1. We also assume that V is instantaneously 
set to zero after firing, i. e., we neglect delays in signal 
transmission between neurons 101. 
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FIG. f. Illustration of the model dynamics. A single 
neuron i undergoes deterministic integrate-and-fire dynam- 
ics (upper right). When this neuron fires, its voltage Vi is 
set to zero, while simultaneously the voltages on all its in- 
hibitory-coupled neighbors are reduced by A (lower right). 



The meaning of the inhibition is illustrated in Fig. ^. 
When a given neuron fires, it instantaneously trans- 
mits an inhibitory action potential to K randomly-chosen 
neighbors whose voltages are each reduced by an amount 
A. This inhibition delays the time until these inhibited 
neurons can reach threshold and ultimately fire them- 
selves. The neighbors of a given neuron are selected at 
random from among all the neurons in the network and 
they are chosen anew every time any neuron fires. Thus 
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the coupling in the network is annealed. While a net- 
work with fixed quenched coupling is biologically much 
more realistic, annealed coupling means that a rate equa- 
tion provides the exact description of the dynamics. For- 
tunately, the annealed and quenched systems appear to 
be statistically identical when the number of inhibitory- 
coupled neurons is sufficiently large. We shall consider 
only this limit in what follows. 

Within the rate equation approach, the distribution of 
neuronal voltages in our inhibitory network is described 
by linear dynamics except at the isolated times when 
a neuron fires. The underlying rate equation admits 
a steady state voltage distribution whose basic proper- 
ties are established analytically in Sec. II. In general, al- 
though the network has a steady response, the dynamics 
of an individual neuron has an interesting history be- 
tween firing events. We study, in particular, the proba- 
bility that a neuron "survives" up to a time t after its 
last firing event in Sec. III. The survival probability de- 
cays exponentially in time with a decay rate that depends 
on the competition between the integration and the in- 
hibitory coupling. This provides a relatively complete 
picture of the dynamics of a single neuron in the net- 
work. A summary is given in Sec. IV and some identities 
are proven in the Appendix. 



II. RATE EQUATIONS AND THE STEADY 
STATE 

For the rate equation description we assume that the 
number of neurons is large and thus a continuum ap- 
proach is appropriate. We define P(V, t)dV as the 
fraction of neurons whose voltage lies in the interval 
(V, V + dV). Then the probability density P(V, t) obeys 
the master equation, 



d_ 

dt 



d_ 

av 



P(V,t)=P 1 (t)6(V) (1) 
+ KP 1 (t) [P(V + A,t)-P(V,t)}, 



where Pi (t) = P(V — l,t). The second term on left-hand 
side accounts for the voltage increase because of the de- 
terministic integration. The first term on the right-hand 
side accounts for the increase of zero- voltage neurons due 
to the firing of other neurons which have reached the 
threshold value V max = 1. The second set of terms ac- 
counts for the change in P(V) due to the processes where 
V + A — > V and V — > V — A (we assume A > since 
inhibitory neurons are being considered). 
In the steady state, Eq. (p|) simplifies to 



dP(V) 
dV 



= Pi S(V) + KPi [P(V + A) - P(V)] . (2) 



Equation <^ may be easily solved by introducing the 
Laplace transform, 



V{s)= ( dVP{V)e sV , 

J — oo 



to give, after some straightforward steps, 
^ ~ K (e- sA - 1) + s/Pi ' 



(3) 



(4) 



The unconventional definition of the Laplace transform 
reflects that fact that the voltage is restricted to lie in 
the range [—00, 1]. 

To solve for the Laplace transform, we first note that 
V(0) — 1 due to normalization. Combining this with 
Eq. (|J), we obtain Pi = (1 + A'A) -1 thus completing the 
solution. The final result is 



V(s) 



1 



A(e~ sA -l) + s(l + AA)' 



(5) 



It may be verified by elementary means that this func- 
tion has a simple pole at s = —A, that is, V(s) — 
A/(s + A) + . . ., where A is the root of 



1) - A(1 + AA) = 0, 



(6) 



and A = (1 - e- A )/[A(l + AA)A - 1]. 

The existence of a simple pole in the Laplace trans- 
form implies that the voltage distribution itself has an 
exponential asymptotic tail as V — > — 00, 



P{V) -> Ae 



\v 



(7) 



The limiting behaviors of the decay constant A may also 
be found from Eq. (^) and give 

x , A -1 In [(AA -1 )] when A ^0, (R s 
* 2(AA 2 )-! when A -> 00. [ ' 
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FIG. 2. Typical simulation results for the steady state neu- 
ronal voltage distribution for a network of 25000 neurons, with 
annealed coupling to K — 50 other neurons and A = 1/50. 
The distribution is shown at 10 time steps. 
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To visualize these results, we have performed Monte 
Carlo simulations of an inhibitory neural network whose 
dynamics is defined by Eq. ([!]). A representative result 
is shown in Fig. initial voltages uniformly distributed 
in [0,1]. Essentially identical results are obtained for 
quenched inter-neuron coupling but all other system pa- 
rameters unchanged. After approximately one time unit, 
the system has reached the steady state shown in the fig- 
ure. Other initial conditions merely have different time 
delays until the steady state is reached. In addition to 
the exponential decay of the distribution for the small- 
est voltages, there are several other noteworthy features. 
First, there is a linear decay of the distribution to its 
limiting value Pi = (1 + AA) _1 as V — > 1 from below; 
this is reflective of the absorbing boundary condition at 
V = 1. There is also a sharp peak at V — 0, correspond- 
ing to the continuous input of reset neurons. For small 
A, almost all voltages lie within the range [0, 1]. This fea- 
ture is reminiscent of the Bak-Sneppen evolution model 
|Tof in which the "fitnesses" of most species lie within a 
finite target range, together with a small population of 
sub-threshold species. 

Finally we note that Eq. ([!]) applies even when the 
number of neighbors for a given neuron is not fixed. In 
this case, K may be interpreted as just the average num- 
ber of interacting neighbors of a given neuron. Therefore 
K can be any positive real number. The voltage off- 
set A can also be heterogeneous, e. g., distributed in a 
finite range with some density p(A). With this general- 
ization, the term P(V + A) in Eq. (^) should be replaced 
by / dAp(A)P(V + A). The resulting equation is still 
solvable by the same Laplace transform technique as in 
the homogeneous network and we now obtain a similar 
solution to that in Eq. (|) but with Pi = (1 + A' (A))" 1 , 
where (A) = f dA p(A) A. 

III. EVOLUTION OF A SINGLE NEURON 

In addition to the steady-state voltage distribution 
P(V), we study the time dependence of an individual 
"tagged" neuron in the steady state. Consider, for ex- 
ample, a neuron which last fired at time Tq which we 
set to for simplicity. The fate of this neuron may be 
generally characterized by its survival probability S(t), 
defined as the probability that this neuron has not yet 
fired again during the time interval (To, To +t), irrespec- 
tive of how many inhibitory inputs it may have received 
(Fig. [}]). A more comprehensive description is provided 
by Qk (t) , the probability that the tagged neuron has not 
yet fired and that it received k inhibitory inputs during 
(0, t). Clearly, S(t) = J2k>o Qk(t)- As we now show, the 
survival probability and Qk(t) exhibit non-trivial first- 
passage characteristics. 

Before discussing the survival probability in detail, let 
us first understand why a tagged neuron must even- 



tually fire and why the survival probability decays ex- 
ponentially at long times. In the absence of interac- 
tions, the voltage of a neuron increases deterministi- 
cally with rate 1. On the other hand, inhibitory spikes, 
each of which reduce the voltage by A, occur stochas- 
tically with rate r = KP X = K/(l + KA). This gives 
rA = KA/(1 + KA) < 1 for the rate at which the volt- 
age decreases due to inhibition. Thus, on average the 
voltage increases at a net rate 1 — rA = (1 + AA) _1 . 
Consequently the voltage of a tagged neuron must even- 
tually reach the threshold and fire. 




FIG. 3. Voltage trajectory of a typical neuron following 
a voltage reset. Because each inhibitory spike reduces the 
voltage by A, a neuron which receives k inputs fires ex- 
actly at time r = 1 + kA. The horizontal tick marks are 
at t = 1, 1 + A, 1 + 2A, and 1 + 3A. 
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FIG. 4. Behavior of the survival probability as a func- 
tion of time. This function is constant within each interval 
[l + (fc — 1)A, 1 + fcA] but exhibits an overall exponential decay 
in time. 

From the above argument, the mean lifetime (t) of a 
neuron between firing events is just the inverse of this 
rate. Thus 

(*) = 1 + KA. (9) 

Notice that the density of neurons which are firing, P 1; 
equals the inverse lifetime, as expected naively. Since the 
evolution of the neuronal voltage is a random-walk-like 
process which is biased towards an absorbing boundary 
in a semi-infinite geometry, the neuron survival probabil- 
ity must decay exponentially with time 111]], 
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S(t) oc e~ t/T when t — > OO. 



(10) 



To understand the asymptotic behavior of S(t), it is 
helpful to consider first the voltage evolution of a sin- 
gle neuron which has experienced 0,1,2,... inhibitory 
spikes. It is important to appreciate that a neuron which 
receives exactly k spikes necessarily fires exactly at time 
1 + kA (Fig. ||). Thus the survival probability is a piece- 
wise constant function which changes discontinuously at 
t k = 1 + kA, k = 0,1,2,... . 

We now determine the value of S(t) at each of these 
plateaux. When t < 1, the tagged neuron has no possi- 
bility of firing and S(i) = Sq = 1. To survive until time 
t = 1 + A, the neuron must receive at least one spike in 
the time interval [0, 1]. Since neurons experience spikes 
at constant rate r = KP\, the probability of the neuron 
receiving no spikes in [0, 1] is e~ r . The survival proba- 
bility for 1 < t < 1 + A thus equals S(t) = Si = 1 - e~ r 
(Fig. |]). Similarly to survive until time 1 + 2 A, the neu- 
ron must receive one spike before t — 1 and a second 
spike before t — 1 + A. By writing the probabilities for 
each of these events, we find, after straightforward steps, 

S' 2 = l-e~ r -re- r(1+A) , 1 + A < t < 1 + 2A. (11) 

This direct approach becomes increasingly unwieldy for 
large times, however, and we now present a more system- 
atic approach. 

To this end, we first solve for Qk(t), the probability 
that the tagged neuron has experienced k inhibitory in- 
puts but has not yet fired. For a constant rate r of in- 
hibitory spikes, the probability that the tagged neuron 
has been spiked exactly k times, with each spike occur- 
ring in the time intervals [tj , tj + dtj 

L<j<fe ! 



equals e rt Y\ 1<r ^,, r dU. Therefore, 



1,2,.. 



, k, 



Q k (t) 



r k e- rt f 



1 



kA-t) J f[ dtj 



(12) 



k > 2, the first spike must occur within [0, 1], while the 
remaining (time-ordered) k — 1 spikes can occur anywhere 
within [t\ ,t]. To evaluate the integral in Eq. (|l2|), we may 
first integrate over ti , . . . , tk ■ This integral is 



dU 



dt 



3 ' ' 



dt k . 



(14) 



The domain of the integral is a simplex of size t — ti 
and the integral is just (t — ti) k ^ 1 /(k — 1)1. Finally, we 
integrate over the region < ti < 1 to obtain 



(rt) k -{rt-r) k 
Qk( - > ~ k\ ■ 



(15) 



This expression actually holds for all k > 0. From 
Eq. (|l5|), we then find that Si = 1 — e~ r in the time 
range 1 < t < 1 + A. 

Generally in the time range l + (m— 1)A < t < 1+mA, 
a tagged neuron which has not fired must have experi- 
enced at least m spikes; therefore Qk = for k < m. 
To determine the non-zero Qfe's - those with k > m - 
note that the first m spikes must obey the constraint of 
Eq. (|l^) , while each of the remaining k — m spikes may 
lie anywhere within the time interval t m and t. The lat- 
ter condition again defines a simplex of size t — t m . This 
gives the contribution (i — t m ) h ~ m / (k — m)\ for the in- 
tegral over these variables. By this reduction, the fc-fold 
integral in Eq. ( |l2| ) collapses to the m-fold integral 



kit) 



r k e -rt 



it tm) 



k-r 



(k — m)\ 



JJdtj. 



(16) 



In particular, for k = to, we may write Q m (<) as 
r m e- rt T m (A), where 



T ro (A) = dti dt 2 ... 

JO Jti Jt m -x 



l + (m-l)A 



dt m . (17) 



Here the step function 8(1 + kA — t) guarantees that the 
voltage of the tagged neuron is below the threshold at 
time t. We must also ensure that the voltage is less than 
one throughout the entire time interval (0,i). The nec- 
essary and sufficient condition for this to occur is that 
the voltage is below threshold at each spike event. This 
determines the integration range in Eq. (|l2T) to be 



tj-! < tj < min[t, 1 + (j - 1)A] 



(13) 



Remarkably, this expression has the simple closed-form 
representation (see Appendix) 



T m (A) = 



(1 + mA)' 



so that 



Qm it) 



r m (1 + mA) m - 1 _ rt 



(18) 



(19) 



for j = 1, . . . , k (we set t = 0). 

We now evaluate Qk (t) successively for each time win- 
dow [1 + (j — 1)A, 1 + jA}. Consider first the range 
1 < t < 1 + A. Here a tagged neuron which has not 
fired must have received at least one spike. Consequently 
Qo{t) = 0. Similarly if a neuron receives a single spike 
and survives until t = 1 + A, the spike must have oc- 
curred in the time range [0,1]. Hence Qi = re~ rt . For 



For large to and also forfc > to, the explicit expressions 
for Qk become quite cumbersome; however, they are not 
needed to determine the asymptotics of the survival prob- 
ability. 

We now use our result for Qk to determine the survival 
probability. For the time range [1 + (to — 1)A, 1+mA], we 
substitute Eq. ( p~6| ) in Sit) = J2k> m Qk(t) an d perform 
the sum over k to give 
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/ ■m 
e-^Udtj. 

3 = 1 



(20) 



This neat expression formally shows that the survival 
probability is constant but m dependent inside the time 
interval [1 + (m — 1)A, 1 + mA]. These properties justify 
the notation in Eq. (p0|). 

The integral on the right-hand side of Eq. (^0|) can be 
simplified by integrating over t m to give 



S n 



S 71 



jm— 1 



T 



-l e 



-r[l + (m-l)A] 



(21) 



This recursion relation allows us to express the survival 
probability in terms of the T„'s with n < m: 



S m — 1 y f n T n e 



■r(l+nA) 



n=0 



Since Soc — 0, we rewrite this as 



-r(l+nA) 



(22) 



which is more convenient for determining asymptotic be- 
havior. 

Let us first use this survival probability to compute 
the average time interval (t) between consecutive firings 
of the same neuron. This is 



<*> = 



S(t) dt 



= 1 + A ^ S m 

m> 1 

= 1 + A n r ™ T n e~ r(1+ " A) . 



(23) 



K 



A -> 0, KA = 0(1) 



appears biologically relevant. In this limit and when the 
time t = 1 + mA is large, the series for S m in Eq. ( p^ ) 
is geometric. Hence, apart from a prefactor, the survival 
probability is given by the first term in this series: 



S m oc r m T m e 



— r(l+mA) 



(25) 



Using Eqs. (|25|), (|18[), and Stirling's formula, we deduce 
that S(t) decays exponentially with time, with the relax- 
ation time in Eq. (nfl) given by 



A 



Pi + ln(l-Pi)' 



(26) 



The different behavior of the two basic time scales, r 
and (t) = 1 /Pi , is characteristic of biased diffusion near 
an absorbing boundary in one dimension [ pT| . Here, the 
mean survival time is simply the distance from the parti- 
cle to the absorbing boundary divided by the mean veloc- 
ity v. In contrast, the survival probability asymptotically 
decays as exp(— v 2 t/D), so that r = D/v 2 , independent 
of initial distance. 

It is instructive to interpret these results for our neural 
network, where a single neuron can be viewed as undergo- 
ing a random walk in voltage, with a step to smaller volt- 
age of magnitude A occurring with probability r dt in a 
time interval dt and a step to larger voltage of magnitude 
dt occurring with probability 1 — r dt. For this random 
walk, the bias velocity is v = 1 — rA = (1 + KA)^ 1 = P\, 
This then reproduces (t) = l/v = 1 + KA = l/Pi. More- 
over, the diffusion coefficient of this random walk is sim- 
ply D oc rA 2 . This random walk description should be 
valid when KA ks 1/Px S> 1 or rA — ► 1, so that a tagged 
neuron experiences many spikes between firing events. 
This then leads to r = D/v 2 cx K 2 A 3 . In this diffusive 
limit of Pi — » 0, the limiting behavior of Eq. (|2^) agrees 
with this expression for t. 



n>l 



In the first line, we use the fact that — S is just the prob- 
ability that the neuron fires at a time t after its previ- 
ous firing, and the last line was derived by employing 
Eq. (22). As discussed previously, the average time be- 
tween firings of the same neuron is (t) = 1 + KA. Equa- 
tion (p3) agrees with this iff the identity 



Y j nT n z n = - 



n>Q 



- rA 



with z — r e rA 



(24) 



holds. This is also verified in Appendix. 

Let us now interpret our results in the context of bio- 
logical applications. Typically, the number of neighbor- 
ing neurons K is large while the spike-induced voltage 
decrement A of a neuron is small, so that the total volt- 
age decrease KA is of order one. In other words, the 
limit 



IV. SUMMARY 

We have determined the dynamical behavior of an 
integrate-and-fire neural network in which there is purely 
inhibitory annealed coupling between neighboring neu- 
rons. The same behavior is also exhibited by a model 
with quenched coupling. Our model should be regarding 
as a "toy" , since so many realistic physiological features 
have been neglected. However, this toy model has the 
advantage of being analytically tractable. We have deter- 
mined both the steady state properties of the network, as 
well as the complete time-dependent behavior of a single 
neuron. The latter gives rise to an appealing first-passage 
problem for the probability for a neuron to survive a time 
t after its last firing. This survival probability is piece- 
wise constant but with an overall exponential decay in 
time. 
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Given the simplicity of the model, it should be pos- 
sible to incorporate some of the more important fea- 
tures of real inhibitory neural networks, such as neurons 
with leaky voltages and finite propagation velocity for in- 
hibitory signals, into the rate equation description. These 
generalizations may provide a tractable starting point to 
investigate more complex dynamical behavior which are 
often the focus of neural network studies, such as large- 
scale oscillations and macroscopic synchronization. 

We thank C. Borgers, S. Epstein, N. Kopell, and S. Ye- 
ung for stimulating discussions. We are also grateful to 
NSF grant No. DMR9978902 for partial financial sup- 
port. 
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APPENDIX A: BASIC IDENTITIES 

We use Eq. ( p2| ) to derive the identity (p"8|). First, we 
note that So = 1. By substituting this into Eq. ( |2^ ) we 
obtain 



E 

n>0 



T n z n 



with z 



-rA 



(Al) 



The requirement that Eq. (Al) holds for arbitrary r and 
A leads to unique set of T n 's. To determine these T„'s, 
let us treat z as a complex variable. Then employing the 
Cauchy formula yields 



1 

2ni 

1 
2tt7 

1 
2tt7 



■ dz 



z n+l 

e r z'(r) 



dr 



[z(r)]«+! 
(1 -rA)e r ( 1+ " A ) 



„n+l 



dr 



(1 + nA) r 



(1 + nA) r 



-A 



(1 + raA)" 
(« - 1)! 



Next we verify Eq. (|2J) by repeating the procedure 
which has just been used to check Eq. flAl| ). As above, 
the quantity nT n may be written in the integral repre- 
sentation 
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